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Delays in biological systems may be used to model events for which the underlying dynamics cannot 
be precisely observed, or to provide abstraction of some behavior of the system resulting more com- 
pact models. In this paper we enrich the stochastic process algebra Bio-PEPA, with the possibility 
of assigning delays to actions, yielding a new non-Markovian process algebra: Bio-PEPAd. This 
is a conservative extension meaning that the original syntax of Bio-PEPA is retained and the delay 
specification which can now be associated with actions may be added to existing Bio-PEPA models. 
The semantics of the firing of the actions with delays is the delay-as-duration approach, earlier pre- 
sented in papers on the stochastic simulation of biological systems with delays. These semantics of 
the algebra are given in the Starting-Terminating style, meaning that the state and the completion of 
an action are observed as two separate events, as required by delays. Furthermore we outline how 
to perform stochastic simulation of Bio-PEPAd systems and how to automatically translate a Bio- 
PEPAd system into a set of Delay Differential Equations, the deterministic framework for modeling 
of biological systems with delays. We end the paper with two example models of biological systems 
with delays to illustrate the approach. 

1 Introduction 

The contribution of computer science in the interdisciplinary field of Systems Biology is to provide lan- 
guages, tools and techniques for the description and analysis of complex biological systems. In particular, 
there exist many formal languages, either based on process algebras or term-rewriting systems, worth 
noting: Bio-PEPA HUlQj]], the stochastic 7T-calculus (23] EH US, Bioambients H3, the fc-calculus US, 
LBS ED, the CLS HUSH, to name but a few. 

Biological systems can often be modeled at different abstraction levels. Specifically, a simple event 
in a model that describes the system at a certain level of detail may correspond to a rather complex 
network of events in a lower level description. The choice of the abstraction level of a model usually 
depends on the knowledge of the system and on the efficiency of the analysis tools to be applied to the 
model. 

Delays can appear in a biological system at any level of abstraction. In particular, there are two good 
reasons for considering delays. Firstly, when there is a network of events whose dynamics cannot be 
precisely observed, or measured in terms of kinetic rates and, secondly, when a complex portion of a 
system is to be abstracted by means of a smaller one. In both cases, a delay may represent the time 
necessary for the underlying network of events to produce some result observable in the higher level 
model. In both cases, the state space of the model with delays will be a reduction of the complete one as 
some parts are abstracted by the delays. 

In mathematics, the modeling of biological systems with delays is mainly based on Delay Differen- 
tial Equations (DDEs), a kind of differential equations, obtained by generalizing Ordinary Differential 
Equations (ODEs), in which the derivative of the unknown function at a certain time is given in terms 
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of the values of the function at previous times. This framework is very general and allows both simple 
(constant) and complex (variable or distributed) forms of delays to be modeled. Practically, DDEs have 
been used to describe biological systems in which events have a non-negligible duration (6l [29J or in 
which a sequence of simple events is abstracted as a single complex event associated with a duration 
EH El. 

It is well-known that the analysis of ODEs can become imprecise due to the approximation intro- 
duced by representing discrete quantities with continuous variables when quantities are close to zero, 
and the same problem can arise in DDEs. Thus techniques for performing stochastic simulation of bio- 
logical systems with delays have been defined. The Delay Stochastic Simulation Algorithms (DSSAs) 
HI El El, often exploiting the Gillespie's Stochastic Simulation Algorithm (SSA) of chemical reactions 
|[T6l . permit computation of a time-trace of the non-Markovian stochastic process underlying a model 
with delays. These algorithms permit different interpretation of delays CD El, in particular, it is possible 
to have a delay '-as-duration approach to the firing of reactions, or a purely delayed one. In the former 
|[T][5]], the reactants are removed at the beginning of a reaction and the products are added at its end, 
namely after the delay plus an exponentially distributed time quantity. In this sense, during the time of 
firing of the reaction, the reactants will not be able to take part in other reactions. However, in HI El the 
need for a different interpretation of delays is discussed via an example of the cell-cycle with delays. 
More precisely, it is shown that, for some biological systems, it is necessary that reactants involved in 
a reaction with delay can have other interactions while waiting for the delay to complete. Indeed, the 
latter interpretation, namely the purely delayed approach, is such that the reactants involved in a reaction 
can have other interactions during the firing of the reaction itself. The reactions are hence scheduled 
and fully performed after the delay and the stochastic time quantity have expired, if the reaction is still 
applicable. 

In this paper, we define a process algebra for the modeling of biological systems with delays. More 
precisely, we use constant delays in the DDEs and, for the DSSAs, we take the delay-as-duration ap- 
proach presented in (TJ. These restrictions are reasonable since they permit us to have a simple algebra 
obtained by extending a well-known one, Bio-PEPA |[T0l[TTi . Also, later versions of this algebra may be 
extended to more complex forms of delay and interpretation of delays. 

Bio-PEPA is a stochastic process algebra for the modeling and the analysis of biochemical networks. 
Bio-PEPA is based on PEPA iflTl . a process algebra originally defined for the performance analysis of 
computer systems, and extends it in order to handle some features of biochemical networks, such as 
stoichiometry and different kinds of kinetic laws. A main feature of Bio-PEPA is the ability to support 
different kinds of analysis. In particular, Bio-PEPA models can be analyzed by performing stochas- 
tic simulation based on the Gillespie's SSA [16] and steady state analysis can be performed on the 
Continuous-Time Markov Chain underlying the semantics of a model. Furthermore, Bio-PEPA models 
can be translated into equivalent deterministic models based on ODEs and, finally, they can be model 
checked using the PRISM 1 181 l27l model checker. The Bio-PEPA modeling paradigm is processes-as- 
species rather than processes-as-molecules, as in the Stochastic %— calculus j23) . This choice, in general, 
permits a smaller state space and hence a model whose analysis is feasible. 

In this paper, we enrich the stochastic process algebra Bio-PEPA with the possibility of assigning 
delays to actions, yielding the definition of a new non-Markovian process algebra: Bio-PEPAd. The new 
algebra is based on the same syntax as Bio-PEPA, hence the definition of Bio-PEPAd systems with delays 
can be easily obtained by adding, to a Bio-PEPA system of the target model, the delay specifications. 
Bio-PEPAd contains two issues to tackle model reduction: the use of the level of concentrations for the 
species, as in Bio-PEPA, and the delays, as a new feature. The semantics of the algebra is given in the 
Starting-Terminating style [9], which allows us to observe the start and the completion of an action as 
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two separate events, as required by delays. Following previous work on Bio-PEPA analysis, we outline 
how to perform stochastic simulation of Bio-PEPAd systems using the DSSAs introduced in [1], and 
how to automatically translate a Bio-PEPAd system in a set of DDEs. 

At the end of the paper we present some examples of biological systems described by Bio-PEPAd. 
In particular, we show the semantics of a toy model in order to clarify the ideas underlying the definition 
of the algebra. Also, we encode in Bio-PEPAd a well-known model of the cell-cycle with delays where 
the passage of cells from different phases of the cell cycle is modeled by a delay. Such a model is then 
translated into a set of DDEs which match with the first definition of the model, appearing in |[28l . We 
end the paper with some discussions about the future work we plan. 

The paper is structured as follows: in Section [2] we recall the definitions of Bio-PEPA that we main- 
tain in the definition of Bio-PEPAd. In Section [3] we separately introduce the syntax and the semantics 
of the language. In Section [4] we present analysis techniques for Bio-PEPAd systems based on DDEs 
and DSSAs. In Section [5] some examples of Bio-PEPAd systems are presented and, finally, in Section [6] 
conclusions and future work are discussed. 

2 Bio-PEPA 

Bio-PEPA iflOlfTTl is a stochastic process algebra, based on PEPA ifTTl . for the modeling and the analysis 
of biochemical networks. The operators of this algebra are designed for easily describing biochemical 
networks. Indeed, features such as stoichiometry of reactions and general kinetic laws can be easily 
described in Bio-PEPA models. Furthermore, as already said in the previous section, the algebra supports 
multiple analysis techniques for the defined models. Stochastic simulations, steady state analysis of the 
CTMC, automatic translation in sets of deterministic ODEs and, finally, model checking analysis can be 
performed on Bio-PEPA models. 

The processes-as-species modeling paradigm of Bio-PEPA permits a smaller state space and, conse- 
quently, a model whose analysis is feasible. A model is described by sequential components representing 
species, and by some model components representing their possible interactions. 

In this section we recall the parts of the definition of Bio-PEPA that we will use to define Bio-PEPA 
with delays. We assume a set of action types srf and we start by recalling the syntax of the processes. 

Definition Bio-PEPA processes are defined by the following grammar: 

S ::= (a,K)opS \ S + S \ C 
P ::= p^p | S(l) 

where op = \. | j | | © | 0, a G stf, «Sf is a set of actions and Z, K G N. We denote with 5? the set of 
all possible species specifications, and we denote with 9* the set of all possible well-formed Bio-PEPA 
processes, as defined in iTlOl . 

The components S and P represent species and their possible interactions, respectively. The element C is 
used to define constant processes. 

Bio-PEPA actions are used to model the events (i.e. the reactions) happening in the biological sys- 
tems we model. The prefix terms in this algebra contain information about the role of the species in the 
actions. In particular, for (a, K)op S we have that (a, k) is the prefix, where a G stf is the action type and 
K is the stoichiometry coefficient of the species in the reaction. The prefix combinator "op " represents 
the role of the species in the reaction. In particular, \. indicates a reactant, t a product, © an activator, 
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an inhibitor and a generic modifier. The actions can appear in a summation term S[+S2, whose 
meaning is the classical "choice" of process algebras. 

In Bio-PEPA a discrete concentration level I is associated with each species. During the simulation 
of a system, the concentration of a species S, denoted by S(l) ranges over {0, . . . ,N$}, where N$ is its 
maximum level of concentr ation statically defined to bound the population size. Also, a fixed step size 
h for all the species is defined. This means that, changing the level concentration of a species by one, 
implies a change in h units of concentration of that species. The granularity, as well as the rate functions, 
are defined in terms of the step size h of the concentration intervals. This choice permits us to deal with 
incomplete information in the exact number of elements, and leads to a reduction of the state space as 
there are less states for each component. 

Bio-PEPA supports multiway synchronization, i.e. synchronization can involve more than two com- 
ponents. This makes it easy to model n-ary reactions, whose modeling in dyadic process algebras is not 
trivial. The term Pi ^Pi denotes cooperation between Pi and P% over the cooperation set Jzf ', which 
determines those activities on which the cooperands are forced to synchronise. For action types not in 
Jz? , the components proceed independently and concurrently with their enabled activities. 

A Bio-PEPA model specification is given in terms of systems, where a system is defined as follows. 

Definition A Bio-PEPA system is a 6-tuple (f,^V,Jf,^,Comp,P} where: 

• y is the set of compartments; 

• jY is the set of quantities describing each species; 

• J?T is the set of parameter definitions; 

• & is the set of functional rate definitions; 

• Comp is the set of sequential component definitions; 

• P is the initial process definition. 

Notice that in Bio-PEPA the kinetic characteristics of the actions are not specified in the syntax 
of processes as in other calculi but, instead, they are separately represented in the notation of system. 
Indeed, in this definition the information about rates is given by & and that about kinetic constants is 
given by Jff, while the initial process definition is P. 

The semantics of Bio-PEPA is given by a Structural Operational Semantics (SOS) [22], similar to 
the one for PEPA. The semantics is based on a capability relation which supports the derivation of 
quantitative information and which is auxiliary to a stochastic relation. The stochastic relation associates 
the rates with the actions performed. The rates are obtained by evaluating the functional rate associated 
with the action, divided by the step size, and by using the quantitative information derived from the 
capability relation, as explained in iTTOl . The use of two relations allows for the association of the rate 
with the last step of the derivation representing a given reaction, which makes it easier to derive the rate 
in the appropriate way, especially in the case of general kinetic laws different from mass-action. 

For the precise definitions and explanations of the components of a Bio-PEPA system, as well as for 
the formal definition of the SOS of Bio-PEPA, we refer to 0jD. 

3 Bio-PEPAd: Bio-PEPA with delays 

In the following sections we separately present the syntax and the semantics of Bio-PEPA with delays 
(Bio-PEPAd). 
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3.1 Syntax and process configurations 

Processes of Bio-PEPAd are defined by the same syntax as Bio-PEPA processes, hence it will be possible 
to easily encode a Bio-PEPA system in one with delays. 

As in Bio-PEPA the general kinetic information is specified separately from the syntax of processes. 
The delays, which are also properties of the actions which can be performed, are similarly represented 
separately in Bio-PEPAd. Indeed, they are defined by functions as 



such that a (a) denotes the delay of action a G stf . From the biological perspective, the choice of using a 
to specify the delays implies that, for every participant in an action a, a unique delay a(a) corresponds, 
which is sound since for each species involved in the reaction modeled by a the delay is unique. A 
Bio-PEPAd system is defined as an extension of a Bio-PEPA one as follows. 

Definition A Bio-PEPAd system is a 7-uple (y,^y,J^,^,Comp,a,P) where: 

• (y, .yy,.Xr, Comp,P) is a Bio-PEPA system; 

• a is a function satisfying £T|) and used to specify the delays of the actions. 
We denote with the set of all possible Bio-PEPAd systems. 

Again, moving from a Bio-PEPA system specification to a Bio-PEPAd one is straightforward. This will 
permit us, in the future, to reuse the system specifications for Bio-PEPA in the context of Bio-PEPAd. In 
order to define the semantics of Bio-PEPAd we define a notion of process configuration. 

Definition Bio-PEPAd process configurations are defined by the following syntax: 



where L is a list of 4-tuples (I, K, a, op) with I, K G N, a G srf and op = 1 1 j | | | 0. We denote with 
^ the set of all well-formed processes configurations. 

The notion of well-formed process configuration is straightforward; any process configuration is 
well-formed if, by removing the list L, its corresponding Bio-PEPA process is well-formed. For clarity, 
in the following we denote a generic process configuration as S(l,L). 

A species 5(7, L) is a species with a discrete level of concentration /, like the species 5(7) in Bio- 
PEPA, but which is currently involved in the actions with delay described by the list L. In particular, if 
the list L contains an entry (/', fc, a, op), this means that there are K levels of concentration of species 5 
involved in a currently running action a which fired when the discrete level of concentration of species 
5 was /', its role in this instance of action is described by op. 

Consequently, L is to be considered as a view of the scheduling list used in the algorithms described 
in HI for simulating stochastic models with delays. More precisely, L is a view of only the scheduled 
events which involve elements of species 5. 

In order to define the semantics of Bio-PEPAd, it is necessary to define some auxiliary functions for 
manipulating the scheduling list L. We denote with Q the domain of all the possible tuples of the form 
(I, K,a,op), and with «Sf^ all the possible lists built over @, hence L G Jz?^. We start by defining, in 



cr:^H>{rel|r>0} 



(1) 



C s ::= (a,K)opC s 
C P ::=Cp£$C P I 



I c s +c s 

C s (l,L) 



c 
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a functional style [20], a function (j> : srf i-> =2?^ i-> *3> to extract the first scheduled events with a given 
action name from the list L as follows: 

(j) (XL = match L with 

|[] ->JL; 

I (l,K,a,op) ::xs — > (l,K,a,op); 
\ x:\xs — > § a xs. 

The function value a L is _L if no entries of action a exist in L (i.e. no actions a are currently running), 
otherwise it is the first entry obtained by a left-to-right recursive scan of L. Notice that we assume the 
syntactic priority of pattern matching. 

Now, we define a function £ : jz/ i— >• h-> Jzf^ such that £ a L is a new list obtained by removing 
the first, if any, occurrence of an action a obtained by a left-to-right recursive scan of L. We define £ as 
follows: 

£ a L = match L with 

in []; 

I (l,K,a,op) :: xs — > xs; 
\ x:\xs — > x :: £ a xs. 

As this is an event list, the ordering of insertion of the tuples determines their ordering for extraction. 
The functions and £, together with the classical append function on lists, namely function @, will be 
used to implement a First-In First-Out (FIFO) policy for insertion and extraction of elements in L. 

Furthermore, as in Bio-PEPA we want to keep the state representation of the models finite by using 
some constraints for the starting of actions. Thus, let us denote the scheduled actions in which the species 
S is involved as a product by n L, where K : h4 5£® is a recursive function defined as 

K L = match L with 



| (I, ?c,a,t) " xs ->• (/, K,a,t) :: nxs\ 
| x :: xs — > % xs. 

The species S(l,L) is currently involved in the delayed actions as follows: for the scheduled actions in 
7t L it is involved as a product, and for the other ones it is involved either as a reactant, a modifier, an 
activator or an inhibitor. Furthermore, let us denote by p : i-> N the function 

p L = match L with 

|[] -► 0; 
(l,K,a,op) ::xs — > K + pxs. 

This function computes how many levels of concentration are involved in all the actions described in its 
input list, regardless of the role of the species in the scheduled event. By following the delay-as-duration 
of approach HI in the interpretation of the delays this implies that, for species S, there are exactly p % L 
levels of concentrations of species S which are currently waiting for their delay to expire before becoming 
available in the species S. These two functions will be used to define the constraints to keep the state 
space finite, as presented in the next sections. 
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A Bio-PEPAd system specification is typically given in terms of a process P G 8? whose semantics 
is given in terms of its equivalent process configuration Pq € . Intuitively, we want the initial term P to 
be modified in the corresponding initial configuration Pq where every species declaration S(lo, []) in Pc 
is such that S(Iq) is in P. The initial process configuration is obtained by adding an empty scheduling list 
to each species because, in the initial configuration, there are no instances of actions with delay currently 
running. Formally, we define, by structural recursion on the processes' syntax a function pt : 8? \- > 
such that 



We augment the definition of Bio-PEPAd systems to 7-tuples of the form (Y, Jf , J^, & ,Comp, G,Pc) 
where Pc is a process configuration of a process. In the following, we may use the notation P to refer 
to either a process or a process configuration; it will be clear from the context to which of them we are 
referring. We denote the extended set of all Bio-PEPAd systems with process configurations as 8? . 

Similarly to Bio-PEPA where the SOS is defined by means of two relations, in this algebra the SOS, 
given in a Starting-Terminating (ST) style 0, is defined by means of three relations that we present in 
the following section. 

3.2 A Structural Operational Semantics 

In the following subsections we define a start relation on process configurations which, in the same style 
as the Bio-PEPA capability one, contains the quantitative information needed to evaluate the functional 
rates and modifies the process configurations to model the start of an action. Also, we define a completion 
relation on process configuration which describes the termination of an action. Finally, along the line 
of the stochastic relation in Bio-PEPA, we define a stochastic relation for Bio-PEPAd systems, based on 
the start and completion relations, which associates rates with transitions. 

The start relation 

This relation contains the quantitative information to compute rates of starting actions. Also, this relation 
modifies the process configuration to model the starting of an action. 

The start relation is — >j,C ^ x @ + x ^ where 6 + G ® + contains, like the capability relation in 
Bio-PEPA, the information to evaluate the functional rate. We define the labels 6 + as 



where w is defined as w ::= [S : op(l, k)} \ w@w, with S 6 y, op an operator, / the level and K the 
stoichiometry coefficient of the components. The list w is of the same type as the one exhibited as a label 
by the capability relation of Bio-PEPA. The Bio-PEPAd start relation is defined as the minimum relation 
satisfying the rules presented in Figure [TJ 

Formally, if a species ((a, k)\S)(1,L) is involved as reactant in an action, then by following the 
delay-as-duration approach fl] its concentration level is decreased by k. Differently, in the case of a 
species involved as a product, its concentration level is not changed because, as previously stated, this 
relation models the starting and not the completion of an action with delay. In the case of a species taking 
part in the reaction as a modifier, an inhibitor or an activator, its concentration level is not changed, as 
expected. Independently of the role of a species, its scheduling list L is modified to record that some 



;U((a, K)op S) = (a, K)op S 

n(s l +s 2 )=s l +s 2 
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({(X,k)IS)(1,L) {a+ ' [S '- U '' K)] \ st {S){1 - K,L@[(l,K,CC,l)}) K<l<N 
((a, K)t S)(l,L) (0t+ ' [5:t( ^ )]) ) rt (5)(Z,L@ [(/, ic, a, t)]) <l + p n L<N 
((a, fc)op 5)(/,L) ( " + ' [5:Qp(/ '' C)]) > rt (S)(Z,L@[(Z, K.a.op)]) 1 < / < A 7 and op = Q,Q 

( a + \s-od!1 k)]) 

((a,K)opS)(l,L) ( - ' - ' ;| > , f (5)(/,L@[(/,K:,a,op)]) k <l <N and op = ® 
Sl Q, L ) ^% st S[(l>,L>) 5 2 (/,L) ^> rt S 2 (/',L') 



(5! + 5 2 )(Z,L) ^>, t 5;(/',L0 (5! +5 2 )(Z,L) ^>, t S 2 (/',Z/) 

A A A' A « A ^ A « A i^lM, p; ^ A 2 

Figure 1 : The start relation — > st C x @ + x ^. 

of its levels of concentration are currently performing action a. Notice that, in order to maintain the 
FIFO property on the scheduling list L, we simply use the append function @ . This is possible because 
of both the multiway synchronization of Bio-PEPAd and the use of fixed deterministic delays. More 
precisely, two instances of the same action starting in two subsequent instants, are assumed to terminate 
in two subsequent instants. This is true in a framework where delays are deterministic however, if they 
were stochastic, the two instances could have multiple orderings for completion. Indeed, because of the 
multiway synchronization in the scheduling list L the two instances will appear subsequently and, hence, 
will complete subsequently. Notice that, in a process calculus with dyadic synchronization, this would 
not have held by simply using function @ . 

We use constraints on the levels to have a finite state space as in Bio-PEPA. The constraints for 
starting the actions are the same as those in Bio-PEPA except the one for the products. In particular, 
the constraints which must be satisfied by a species S(l,L) to fire an action as a product is, as expected, 
< / + p 7t L < N, if ./V is its maximum level. Intuitively, this means that the levels of concentration in the 
state, /, plus those which already scheduled to be produced, p n L, must not cross the capacity threshold 
N. 

The starting of the action a, in the style of the ST semantics J9], is denoted by the action symbol 
a + , exhibited as a label for all the start derivations. The composition of the derivations of this relation is 
straightforward. 

Some further considerations and comparisons with Bio-PEPA are useful. Firstly, when the actions 
have no delay as in Bio-PEPA, whenever an action fires, the changes in the process are immediately 
visible in a one-step derivation, since the Bio-PEPA capability relation modifies the process according 
to the action. In this algebra, as the instant in which an action starts and terminates are detached, then 
the start relation modifies the process to represent just the starting of the action. Indeed, another relation, 
which does not exist in the semantics of Bio-PEPA, will model the termination of a currently running 
action. 

Secondly, by comparing the algorithm presented in [ 1 ] and the definition of this relation, it is clear 
that the modification of the process to reflect the starting of an action corresponds to scheduling of the 
reaction in the scheduling list. 
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a L = (/, fc,a,t) 



(j) ctL= (l,K,a,op) 



S(l>,L) ^m co S(V +k ,C0CL) S(l>,L) ^^A co 5(/ , C „ L) 



0/7 = 1,0,0, e 



5 2 (/,L) i^l> C0 5 2 (/',L' 



(Si +S 2 )(l,L) ^% co S[(V,L>) (S 1 +S 2 )(l,L) ^% co S> 2 (V,L>) 



*1 2 'Co *i 



A g P 2 " ! ^rn Pll W A 



co 1 J^f 2 



Figure 2: The completion relation — > C . C ^x8 x ^ . 
The completion relation 

This relation is used to model the completion of an action with delay which is currently running. Also, 
this relation contains quantitative information needed to re-compute the functional rate of the action at 
the moment in which it started. 

The completion relation is — y co C ^ x x ^ where 6 is a similar label to the one exhibited 
by the start relation. We define the labels 6 as 

6 := (a~,w) 

where w is defined as for the start relation. The completion relation is defined as the minimum relation 
satisfying the rules of Figure [2] 

Formally, for a species S(l,L) it is possible to get the instance of a currently running action a, if any, 
by applying function (j). More precisely, this permits us to get, from all the possible instances of actions 
a, the first which has been scheduled, <p a L, and, hence, the first which will terminate. If the species 
is involved as a product, then it is necessary to increase, as defined by the delay-as-duration approach, 
its concentration level by adding the scheduled products. Otherwise, whatever the role of the species, 
its concentration level must remain constant. Independently of the role of the species in the action, the 
scheduling list is modified by means of the function £, hence a new list £ a L is produced by removing 
from L the entry which was computed by function . 

Obviously, no constraints are stated for the completion of a currently running action, as the appropri- 
ate ones are checked before the starting of actions. 

The completion of the action a, in the style of the ST semantics (9), is denoted by the action symbol 
a~, exhibited as a label for all the completion derivations. The other label, namely the list w, is defined 
like the one exhibited by the start relation. The composition of this relation with the other operators is 
straightforward and very similar to the composition of the derivations of the start relation. 

Some further considerations are worth noting. Firstly, this relation is a new one with respect to 
the Bio-PEPA semantics. Again, in a framework where actions have no delays the contribution of this 
relation to the semantics would have been given by means of a unique relation. Also, as the role of this 
relation is to model the completion of an action, it chooses actions to terminate from the list which is 
associated with the species, namely the list of actions currently running. The start relation, differently, 
chooses the action to fire from the species definition. 

Furthermore, as we want the completion relation to exhibit quantitative information to recompute the 
functional rate of the action at the moment at which it started, then the labels exhibited by this relation 
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{f,Jf,J(f,&,Comp,a,P 




> a {f \JV ', Jf \& \Comp,a,P' 




{y,^,J^,^,Comp,a,P 



(a- ,r a [w,jV \X],a(a)) 



> a {V,jY,J(f,&,Comp,o,P' 



Figure 3: The stochastic relation — > S C & x T x 



are very similar to those exhibited by the start relation, even if they are computed starting from a L. 
This permits us to have a unique policy for computing the functional rates from the input lists, obtained 
by derivations of the transitions of these relations. 

The stochastic relation 

The stochastic relation permits us to associate rates to transitions. Also, this transition permits us to 
observe changes in a Bio-PEPAd system due to either the stalling or the completion of an action. 

The stochastic relation is -K s C^xfx^ 1 where y € T is defined as 



with a € s/, r a € R + and o a € R. As in Bio-PEPA, r a represents the parameter of an exponential 
distribution and, as expected, all activities enabled attempt to proceed but only the fastest succeeds. 

As this relation is defined on the set namely the set of all possible Bio-PEPAd systems with pro- 
cess configurations, whenever we refer to the semantics of a system (V \JV , \& ,Comp,a ,P), where 
Pis a process, we assume we apply the stochastic relation to the system {"f \JY , \Comp, a,/i(P)). 
Again, this is necessary because P is not a process configuration, and we want to build, from P, the 
corresponding initial configuration ju(P), and then we want apply the semantics to the system. 

The stochastic relation is defined as the minimum relation satisfying the rules given in Figure [3] 

Formally, the starting of an action a, obtained by composition with a derivation of the start relation, 
is denoted by symbol a + . The completion of an action is obtained by composition with a derivation of 
the completion relation, as denoted by symbol a . 

The rate of any action is computed as in Bio-PEPA, namely as r a [w,^V,J^] = f a [w,^V,J^] ■ h~ l . 
For the explanation of how the rates are computed because of the levels we refer to ifTOj . For any possible 
derivation of the stochastic relation, the value a {a) denotes the delay of the action a. 

A Stochastic Labelled Transition Systems can be defined for a Bio-PEPA system with delays. 

Definition The Stochastic Labelled Transition Systems (SLTS) for a Bio-PEPAd system is (^,r,— y s ) 
where — » v is the minimal relation satisfying the rules given in Figure [3] 

4 Analysis techniques 

In this section we present some analysis techniques for Bio-PEPAd systems along the line of those 
presented in |[T0l for Bio-PEPA systems. Firstly, we discuss the automatic translation of a Bio-PEPAd 
system into a set of Delay Differential Equations (DDEs). Secondly, we discuss how to apply a Delay 



y:= (a + ,r«,a a ) | (a ,r, 
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Stochastic Simulation Algorithm (DSSA) to compute the stochastic time-evolution of a Bio-PEPAd 
model. 



4.1 Translation in Delay Differential Equations 

Whenever phenomena presenting a delayed effect are described by differential equations, we move from 
ODEs to DDEs. In DDEs the derivatives at current time depend on some past states of the system. The 
simplest form of DDE considers constant delays, namely consists of equations of the form 

^ = f x (t,X(t),X(t - ffi), • • • ,X(t - O n )) 

with d > . . . > G n > 0, Gi G M and X(t — a,) denotes the state of the system at the past time t — a,. 
This form of DDE allows models to describe events which have a fixed duration. Hence it is natural, in 
the context of Bio-PEPAd, to reason about the translation of a model into a set of DDEs. Furthermore, 
similar work has been presented in [10] for translating a Bio-PEPA system into a set of ODEs. 

In order to define the encoding it is important to recall that we defined Bio-PEPAd in terms of Bio- 
PEPA. This means that, given a system specification ( 'f , jV , , , Com p,o,P) where P is a valid Bio- 
PEPA process, we just need to modify the algorithm defined in [10] to add the information provided by a 
concerning the delays. Formally, the results for Bio-PEPA permit us to encode \JY , ,Comp,P) 
in a set of ODEs by using the definition of the stoichiometry matrix associated with P. 

The algorithm presented in iflOl consists of three steps. In the first the stoichiometry matrix is defined, 
in the second the kinetic law vector Vkl is derived and in step three the deterministic variables are 
associated with the components. Steps (1) and (3) are unaffected by the use of delays; hence we preserve 
them. 

Step (2), namely the definition of the kinetic law vector Vkl, must be changed. Such a vector contains 
the kinetic law of each reaction; we will explain the definition via an example. For instance, for an action 
a involving species Si and 52, and with mass action kinetics /ma(^), its original entry in the vector Vkl 
would be kxs { {t)xs 2 {t) where xs l and xs 2 are the deterministic variables representing species Si and S2 
species, respectively. The variables depend on the state at time t, but in the context of DDEs, the delays 
of the actions become dependencies on the past states of the system. Hence, for that particular example, 
the correct entry in the vector v KL must be kx S{ (t — o(a))xs 2 (t — o(a)). Step (2) can be generalized 
adding, in the process of the definition of Vk l , the delays of the form x$(t — cr(ce)) for all species S and 
actions a. 

The DDE system can be defined in the same way as the ODE one, namely as dx/dt = D x Vkl where 
x and D are the results of step (3) and (1) of the algorithm, respectively. The initial conditions are, 
however, different from the ones defined for ODEs. In particular, the DDEs, because of the delays, must 
be defined also in the interval [to — a(a);?o] where a is the action with maximum delay. 

It is not possible to define a universal initial condition for the DDEs systems as every possible config- 
uration will affect the dynamics of the whole system. Sometimes the initial conditions of a species S are 
defined via a constant function 0s(7) for t € [to — a (a); to] such that <j>s(t) = hls,o where ls$ is the initial 
concentration level for S in the Bio-PEPAd model and h is the step size for the concentration levels (see 
[1J for details). In general, we leave this part of the translation to the modeler who will tune the initial 
conditions with respect to the specification of the target system. 



96 



Modeling biological systems with delays in Bio-PEPA 



4.2 Stochastic Simulation 

The stochastic simulation of biological systems is typically based on the SSA by Gillespie lfl6l and its 
variants. Anyway, the SSA, as well as all its variants, are not able to deal with actions with delays but 
only with Markovian actions. As a consequence, some DSSAs UlEl have been defined to perform the 
stochastic simulation of a system where actions have a fixed delay. 

In this section we briefly explain how to perform the stochastic simulation of a Bio-PEPAd system by 
using the DSSAs presented in (H, where all the reactions follow a delay-as-duration approach. Also in 
this context the choice of reusing part of the Bio-PEPA definitions in Bio-PEPAd is crucial. In particular, 
this permits us to completely re-use the techniques defined in iflOl to perform the stochastic simulation 
of Bio-PEPA systems by using the SSA. 

The main steps in preparing a Bio-PEPAd system for the application of the DSSA are the two. 
Firstly, given an initial system (Y, JV , .Jtf , ^,Comp, G,P), where P is a Bio-PEPA process as well as a 
Bio-PEPAd one, a vector describing the initial numbers of molecules to be simulated must be obtained 
by an encoding of P. Secondly, the actual rates of the reactions have to be defined, by cases, starting 
from 

Both the two issues are independent of the target algorithm (hence of the delays). More precisely, 
since both the SSA and the DSSA simply assume a vector X(/o) describing the initial state of the system, 
and as P is also a Bio-PEPA process since it is not in a configuration, then we can simply use the 
techniques developed in ifTOll to compute X(?o) from P. 

Also the definition of the actual rates of the reactions can be done again in the same way for both 
Bio-PEPA and Bio-PEPAd systems. This is possible since the DSSA we refer to is based on the SSA. 
Indeed, it assumes the same type of definition of propensity functions to compute the probability of 
the reactions,. Finally, as is defined in the same way for both Bio-PEPA and Bio-PEPAd, then the 
techniques developed in IfTOll to define from J£" the actual rates can be used also in the context of Bio- 
PEPAd. 

Once these two steps have been performed, the resulting systems can be simulated by the DSSA 
presented in 0] |H where all the actions follow a delay-as-duration approach. 



5 Examples 

In this section we present a toy example of a Bio-PEPAd model to illustrate the approach, followed by 
an encoding of a well-known model of the cell cycle in Bio-PEPAd. 

A toy example 

In order to clarify the modeling with Bio-PEPAd we present a toy example of a model. We assume a 

k a' 

single reaction channel of the form A — — > B to denote a transformation of an element of species A into 
an element of species B with a kinetic constant k and a delay a'. The initial state contains three elements 
of species A and no elements of species B; formally it is described by the vector X(? ) = (3,0) r . 
The Bio-PEPAd processes modeling the species can be easily defined as follows: 

A d = (a, l)]A B=f{a,\)tB 

where a is the action which models the reaction, the functional rates are defined according to the mass- 
action kinetics, f a = /ma(^) and the delay is defined according to the function a(a) = a'. The Bio- 
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Figure 4: A graphical representation of the stochastic process obtained by applying the semantics to the 
process configuration for the toy example. 

PEPAd process describing the interacting components is A 9*\B. By considering levels we assume the 

jo;} 

species to have some maximum levels Na and Nb where Na > 3 and Nb > 3. The initial levels of 
concentrations are described by the vector X(?o), and the initial configuration of the process, obtained by 
applying function pL, is the following 

A(3,[])g}B(0,[]). 

By applying the stochastic relation to the system with this process configuration we obtain all the 
possible evolutions of the configuration. The obtained LTS, as expected, is finite, and, because of the 
delays, it corresponds to a non-Markovian stochastic process. Intuitively, there is a one-to-one corre- 
spondence between both the states and the transitions of the LTS and of the stochastic process, exactly 
as between the LTS for Bio-PEPA systems and the CTMCs. 

A graphical representation of the stochastic process is given in Figure [4] In that figure, all the states 
are represented as circles where the notation (n\,n2):m represents the discrete levels of concentration 
of the species A, n\, and B, ni. The number m represents the number of instances of the unique possible 
action a currently scheduled in the state. All the arrows represent stochastic derivations of the whole 
system, where the labels are exactly those computed by that relation. The full arrows represent stochastic 
derivations based on start derivation, empty arrows represent stochastic derivations based on completion 
derivation. For this particular example, any empty arrow built from a derivation with a rate r refers to the 
completion of the unique action started with the same rate r. 

Figure [5] presents a table showing the explicit mapping of the states described in Figure |4] and the 
corresponding process configuration obtained by the semantics. For the sake of clarity, as in this simple 
example there is just one action, a, and A always participates in that action as a reactant and 6 as a 
product, this information is omitted from the scheduling lists. 

As expected, this system, starting from the initial configuration X(?o), namely state(3,0) : 0, eventu- 
ally reaches the final state (0, 3) : 0, which corresponds to the final configuration A (0, [ ]) r^B(3, [ ]) and 
to the vector X(f') = (0,3) r , for some t' > t . 



A model of the cell cycle with delays 

In this section we encode in Bio-PEPAd a model of the cell cycle with delays as presented in [lj. Such 
a model is obtained by simplifying a DDEs model of tumor growth that includes the immune system 
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state 


process configuration 


(3,0) :0 


A(3,[])g}2?(0,[]) 


(2,0) : 1 


A(3,[(3,1)]){*B(0,[(0,1)]) 


(2,1) :0 


A(2,[]){J } B(1,[]) 


(1,0) :2 


A(l,[(3,l),(2,l)])g|B(0,[(0,l),(0,l)]) 


(1,1): 1 


A(l,[(2,l)])g}B(l,[(l,l)]) 


(1,2) :0 


A(l,[]){*2?(2,[]) 


(0,0) : 3 


A(l,[(3,l),(2,l),(l,l)])g}B(0 ) [(0,l),(0,l),(0,l)]) 


(0,1) : 2 


A(l,[(2,l),(l,l)])g|B(l,[(l,l),(l,l)]) 


(0,2) : 1 


A(0,[(l,l)])g}fi(2,[(2,l)]) 


(0,3) :0 


A(0,[]){J } B(3,[]) 



Figure 5: A table stating the correspondence between the states represented in Figure |4] and the process 
configurations obtained by the semantics. 



response and a phase-specific drug able to alter the natural course of action of the cell cycle of the tumor 
cells EH . 

The model of the cell cycle with delays has been analyzed in [ Q in order to discuss two possible 
interpretations of delays in the delay stochastic simulation algorithms, a delay-as-duration approach and 
a purely delayed approach. In this section, we simply show how to encode that model in Bio-PEPAd and 
for a detailed analysis of the model we refer to that paper. 

The cell cycle is a series of sequential events leading to cell replication via cell division. It consists 
of four phases: Gi, S, G2 and M. The first three phases (Gi, S, G2) are called interphase. In these phases, 
the main event which happens is the replication of DNA. In the last phase (M), called mitosis, the cell 
segregates the duplicated sets of chromosomes between daughter cells and then divides to form two new 
cells in their interphase. The duration of the cell cycle depends on the type of cell (e.g a normal human 
cell takes approximately 24 hours to perform a cycle). Cell death via apoptosis may happen in any phase 
of the cell cycle. 

The Bio-PEPAd model considers two populations of cells: 7}, the population of tumor cells during 
cell cycle interphase, and Tm, the population of tumor cells during mitosis. We consider four possible 
actions, a, /3, 7 and 8, one for each of the events that we want to model. In particular, action a models 
the passage from the interphase to the mitotic phase, with rate a\, j8 models the mitosis, with rate a\, 7 
the death of a cell in the interphase, with rate d2, and 8 the death of a cell in the mitotic phase, with rate 
d^. All the rates in the model refer to mass action kinetics. 

The Bio-PEPAd model is defined by the following species definitions: 

r / t / (a,l)|+(j3,2)t+(7,lH 
T M = (0,1)1+03, 1H + 



where the species behave as reactants or products, depending on their role as previously specified. Also, 
as all the actions obey a mass action kinetic law, we simply assume f a = /ma(oi), fp = /ma (04), fy = 
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/ma (d 2 ) and fg = f MA (rf 3 ). The Bio-PEPAd process modeling the interactions is given by 

T M { ^ } T M (n M ) 

where Rq and represent the initial concentration levels for the cells in the interphase and in the mitotic 
phase, respectively. Notice that y and 8 are not in the cooperation set since model reactions involving a 
single species. Also, we recall that this is also a valid Bio-PEPA process specification. 

A delay a' is used to model the duration of the interphase, hence the passage of a tumor cell from 
the population of those in the interphase to the population of those in the mitotic phase, namely the event 
modeled by action a, is delayed. To specify the delay in the Bio-PEPAd system to analyze, it is enough 
to define a function a where 

a(a) = a' o(fi) = a(y) = a(<5) = 0. 

As a consequence, the Bio-PEPAd process initialized by applying function pt, namely the process con- 
figuration 7}(?1q, [ ]) r^s-, Tm{tiq , [ }), together with the function a, completes the definition of the Bio- 
PEPAd system representing the cell cycle model. 

By applying one of the techniques discussed in this paper this system can be analyzed. In particular, 
the Bio-PEPAd model can be automatically translated into a set of DDEs by applying the algorithm 
presented in Section |4~T1 By computing the following vector of the kinetic laws 

v KL = {a x Ti(t - cr(a)), a 4 T M (t - a(/3)), d 2 T,{t - <r(y)), d 3 T M (t - a(8))) T 
= {aiT^t-o'), a 4 T M (t), d 2 Ti{t), d 3 T M (t)) T 

the following set of DDEs can be computed: 

— - = laa^M — d 2 Tj — a[T/(t — o') — — = a\Ti(t — o') — d^M — cl\Tm ■ 
at at 

As expected, this DDEs system is analogous to the one presented in (H. The terms d 2 Tj and d^M 
represent cell deaths. The cells reside in the interphase at least a' units of time; then the number of cells 
that enter mitosis at time t depends on the number of cells that entered the interphase a' units of time 
before. This is modeled by the terms Tj(t — a') in the DDEs. Also, each cell leaving the mitotic phase 
produces two new cells in the 7} population, as given by terms —a^M and 2a 4 TM- As a consequence, 
by defining the appropriate initial conditions for the resulting DDEs system it would be possible to 
reproduce the results presented in HI for the deterministic model. 

As far as the stochastic analysis of the Bio-PEPAd systems is concerned, we can notice that the 
system we defined corresponds to the following set of reactions 

7} -^4 Tm with delay a Tm -4 27} 7} -4 Tm ■ 

Again, this is exactly the same reactions-based model used in []T| to compare the deterministic and the 
stochastic models for the cell cycle. Consequently, by applying the DSSA as explained in Section l4~2l it 
would be possible to reproduce the results presented in [l] for the stochastic model. 
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6 Discussion and conclusions 

In this paper, we have enriched the stochastic process algebra Bio-PEPA with the possibility of assigning 
delays to actions, yielding the definition of a new non-Markovian process algebra: Bio-PEPAd. The use 
of delays in biological systems is suitable to model events for which the underlying dynamics cannot be 
precisely observed. Also, delays can be used to abstract portions of systems, leading to a reduced state 
space for models. From this point of view Bio-PEPA, which is based on the idea of levels to tackle the 
problem of state space explosion, was an appropriate candidate for defining our algebra. 

The algebra is based on the syntax of Bio-PEPA. Hence the definition of Bio-PEPAd systems with 
delays can be easily obtained by adding, to a Bio-PEPA system of the target model, the delay specifica- 
tions. 

The semantics of the firing for the actions with delays is the delay-as-duration approach, as presented 
in the definition of DSSAs. In future work, we may enrich Bio-PEPAd with the other interpretation of 
delays presented in jTlEl, in order to have the purely delayed approach and its combination with the one 
we currently consider. 

The semantics of the algebra has been given in the Starting-Terminating style. This permits us to 
observe the start and the completion of an action as two separate events, as required by delays. In future 
work, we will consider equivalence relations for Bio-PEPAd systems and processes, as done in lfl4l for 
the Bio-PEPA ones. 

In keeping with the techniques developed for analyzing Bio-PEPA models, we outlined how to per- 
form stochastic simulation of Bio-PEPAd systems and how to automatically translate a Bio-PEPAd sys- 
tem in a set of Delay Differential Equations, the deterministic framework for the modeling of biological 
systems with delays. Moreover, the software framework for Bio-PEPA [7 ] could be extended to provide 
a tool for the automatic analysis of Bio-PEPAd systems. 

As a proof of concept, we presented two examples of Bio-PEPAd systems. The first one, a toy 
example, has been shown to illustrate the semantics we defined. The second one, a well-known model 
of the cell-cycle where phase passages are abstracted by means of a delay, has been presented in order to 
show the translation of the Bio-PEPAd system into both a stochastic process with delays to be simulated 
by a DSS A, and a set of DDEs which can automatically derived by the system specification. 

In the future, we plan to define Bio-PEPAd models of biological systems with delays and to analyze 
such models using the anlysis techniques we defined in this paper. 

An interesting area for further future work will be to compare Bio-PEPAd with non-Markovian 
Stochastic Petri Nets such as DSPN fl5l 
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